Jamming I: A volume function for jammed matter 
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We introduce a "Hainiltonian"-like function, called the volume function, indispens- 
able to describe the ensemble of jammed matter such as granular materials and emul- 
sions from a geometrical point of view. The volume function represents the available 
volume of each particle in the jammed systems. At the microscopic level, we show that 
the volume function is the Voronoi volume associated to each particle and in turn we 
provide an analytical formula for the Voronoi volume in terms of the contact network, 
valid for any dimension. We then develop a statistical theory for the probability dis- 
tribution of the volumes in 3d to calculate an average volume function coarse-grained 
at a mesoscopic level. The salient result is the discovery of a mesoscopic volume func- 
tion inversely proportional to the coordination number. Our analysis is the first step 
toward the calculation of macroscopic observables and equations of state using the 
statistical mechanics of jammed matter, when supplemented by the condition of me- 
chanical equilibrium of jamming that properly defines jammed matter at the ensemble 
level. 
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I. INTRODUCTION 

The development of a statistical mechanics for gran- 
ular matter and other jammed materials presents many 
difficulties. First, the macroscopic size of the constitu- 
tive particles forbids equilibrium thermalization of the 
system. Second, the fact that energy is constantly dis- 
sipated via frictional interparticle forces further renders 
the problem outside the realm of equilibrium statistical 
mechanics. In the absence of energy conservation laws, a 
new statistical approach is needed in order to describe the 
system properties. Along this line of research, Edwards 
[ll has proposed to replace the energy by the volume as 
the conservative quantity of the system. Then a canon- 
ical partition function of jammed states can be defined 
and a statistical mechanical analysis ensues. 

While it is always possible to measure the total vol- 
ume of the system, it is unclear how to treat the volume 
fluctuations at the microscopic level. Thus, it is still an 
open problem the definition of an analogous "Hamilto- 
nian" that describes the microstates of jammed matter. 
Such a "Hamiltonian" is called the volume function 
denoted W. The idea is to partition the granular ma- 
terial into N elements and associate an additive volume 
function to them, W,;, such that the total system volume, 
W, is 
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(1) 



From a theoretical perspective, initial attempts to de- 
fine the volume function involved modelling under mean- 
field approximations proposed by Edwards not taking 
into account the contact network. The necessary defini- 
tion in terms of the internal degrees of freedom has been 
pursued by Ball and Blumenfeld [5, 's'l who have shown by 
an exact triangulation method that the volume defining 
each grain can be given in terms of the contact points us- 
ing vectors constructed from them. The method consists 



of defining shortest loops of grains in contact with one 
another, thus defining the void space around a central 
grain. 

A simpler version for the volume function was also 
given by Edwards as the area in 2d or volume in 3d 
encompassing the first coordination shell of the grains 
in contact [4]. The resulting volume is the antisymmet- 
ric part of the fabric tensor, the significance of which is 
its appearance in the calculation of stress transmission 
through granular packings This definition is only an 
approximation of the space available to each grain since 
there is an overlap of volumes for grains belonging to the 
same coordination shell. Thus, it overestimates the total 
volume of the system: ^ Wi > W [4] . 

Furthermore, both definitions of Wi in 0, [1, |3| are 
proportional to the coordination number of the grain. 
This is in contrast to expectation since the free volume 
available to a grain should decrease as the number of con- 
tacts increase. Indeed, this observation is corroborated 
by experimental studies of jammed granular matter using 
X-ray tomography Q to determine the volume per grain 
versus coordination number. 

Based on the idea that the volume function represents 
a free volume available per grain, we introduce a new 
"Hamiltonian" to describe the microstates of jammed 
matter. We analytically calculate the volume function 
and demonstrate that it is equal to the Voronoi volume 
associated to each particle, partitioning the space into 
a set of regions, associating all grain centroids in each 
region to the closest grain centroid. Even though the 
Voronoi construction successfully tiles the system, its 
drawback in its use as a volume function was that, so 
far, there was no analytical formula to calculate it. Our 
approach provides this formula in terms of the contact 
network. Furthermore we introduce a theory of volume 
fiuctuations to calculate a coarse-grained average volume 
function defined at the mesoscopic level that reduces the 
degrees of freedom to only the coordination number z. 
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We find that the volume function is inversely propor- 
tional to z in agreement with the tomography experi- 
ments of [1]. Our analysis also provides an equation of 
state, relating volume with coordination number in the 
limit of fully random system. Indeed, it predicts with 
good accuracy the limiting cases of random loose and 
random close packing fractions. Our results allow con- 
struction of a statistical partition function from which 
macroscopic observables can be calculated. This case is 
treated in more detail in the second part of this work: 
Jamming II Q. 




Outline 



This paper is the first installment in a series of papers 
devoted to different aspects of jammed matter and is the 
main theoretical contribution for the subsequent statis- 
tical mechanics approach. The present paper is an ex- 
tended version of the Supplementary Information Section 
I in [6*1 . The outline is as follows: Section |TT] details the 
development of the microscopic volume function in terms 
of the particle coordinates. The relation of this form with 
the Voronoi volume of each particle is discussed in Sec- 
tion mil Section IIVI discusses the statistical theory to 
calculate the probability distribution of the Voronoi vol- 
umes leading to the average mesoscopic volume function 
discussed in Section |Vl Section IVll tests the assumptions 
of the theory and we finish with the outlook in Section 

Em 



FIG. 1: Definition of the volume function. Particle i moves 
in the direction of s as indicated such that the deformation 
energy is below a threshold e. When averaged over all direc- 
tions s, this process defines the available volume under the 
threshold e demarcated by the thin dotted line. The volume 
function is then obtained in the limit of infinite rigidity of the 
particles, a. — > oo. 



energy threshold e is introduced in order to drive the 
particle in a certain direction s as indicated in Fig. [T] 
A small displacement br = 6rs of particle i along the 
s direction results in an increase of energy E between 
particles i and its neighbors: 



II. MICROSCOPIC VOLUME FUNCTION 

We start by defining a volume function for rigid spher- 
ical grains of equal size in terms of particle positions. 
The volume function represents the available volume to 
the particle with the constraint of fixed total system vol- 
ume. Since we are dealing with rigid jammed particles, 
the available free volume is in principle zero since the 
particle by definition cannot move. However, we allow 
the particle to move by introducing a soft interparticle 
potential and then taking the limit of particle rigidity 
to infinity. The resulting volume is well-defined, repre- 
senting the free volume associated with each grain in the 
jammed packing. 



(5r ■ r i > 



s-f ij>0 

where the sum is taken over all the neighbors of particle i. 
The soft pair potential fa can be any repulsive function 
provided it approaches the hard sphere limit when the 
control parameter a oo, implying 



lim ^^Oyx<y. 
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(3) 



[A possible function is simply fa{x) = x", the case a = 
5/2 corresponds to the Hertz potential]. This condition 
implies that 



We consider a rigid particle of radius R jammed at po- 
sition ri in contact with another particle at position rj 
such that fij — fj ~ fi — rijfij (see Fig. [T]). In order 
to calculate the volume associated with such a particle 
we allow it to move by introducing a generic interparti- 
cle soft-potential, /q,, determined by an exponent a gov- 
erning the rigidity of the particles (see below). A small 



lim Yfa{x.i) =/oo(max(x.O) x lim E j- / " — ~i — = 
Q^oo I Q-»oo /oQ(maxi(xi)j 

=/oo(max(a;i)), 

i 

(4) 

and therefore from Eq. ([2]), we obtain: 

= /oo(max(s • fy )(5r). (5) 



We define the available volume to a grain, Wf , under 
the energy threshold e as: 



J Q{e^E)dV^ j J Q{e- E)5r'^-^d[6r]d[s], 

(6) 

where d is the dimension of the system. 

The integration over Sr can be simplified when a — > oo 
since for a fixed direction s we have: 

/■OC 

lim / Q{e- E)6r'^-^d[5r]^ 

/>oo 

= hm / Q[e — fa{ max {s-fij)Sr)]Sr'^^^d[Sr] (X 



S-f ij>0 



(X min (s • fijY 
s-fij >o 



Thus, we obtain: 



Wf oc ® min [s ■ rij) ''ds. 



(7) 



(8) 



Equation ([8]) can be interpreted as follows: for each 
direction J, the available volume is determined by the 
particle position whose projection of the distance to par- 
ticle i in the s direction is minimal. The total volume is 
then the average over all directions s. The proportional- 
ity constant in Eq. ([8]) can be determined because Wf 
is equal to the volume of the grains, Vg, when the coor- 
dination number z oo, suggesting that, in this limit, 
mmg.f.^yols ■ fij)^'^ = 1 for any s. That is. 



Wf = 



<f ds 



/min [s-fij) '^ds. 
s-fij>Q 



For mono-disperse spherical particles, Vg 



Thus, we have 
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(9) 
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s-fij>o s ■ r,- 



ds 



s-fij>o 2s ■ fi 



(10) 



ds. 



where we have replaced = 2R for nearest neighbors in 
the last equation. This allows us to generalize the volume 
formula to satisfy additivity and relate it to the Voronoi 
volume as shown below. 

Equation pi7|) is not additive. However, the formula 
becomes additive when considering all particles rather 
than the nearest neighbors in the calculation of the min- 
imum integrand in Eq. (jlOp . This approach is justified 
since non-contact particles may contribute to the energy 
of deformation in Eq. Further, if the minimum in 

Eq. PU)) is taken over all the particles in the packing, Wf 
is exactly equal to the Voronoi volume, which is obviously 
additive. In turn we provide a formula for the calculation 
of the Voronoi volume in terms of particle positions, as 
shown next. 




FIG. 2: The limit of the Voronoi cell (grey area) of par- 
ticle i in the direction s is li{s) — min j lij{s) /2, where 
lij ~ Tijl cos 9ij. Then the Voronoi volume is proportional 
to the integration of li{s) over s as in Eq. (|14|l . 



III. A FORMULA FOR THE VORONOI 
VOLUME 

First, we recall the definition of a Voronoi cell as a 
convex polygon whose interior consists of all points which 
are closer to a given particle than to any other. 

Formally, the volume of the Voronoi cell of particle i 
can be calculated as (see Fig. [2]): 



W7 



drds 



kisfds 



(11) 



where li{s) is the distance from particle i to the bound- 
ary of its Voronoi cell in the s direction. Note that this 
definition is valid whether the particle j is in contact with 
i or not. If we denote the distance from particle i to any 
particle j at s direction as 



(12) 



then li{s) is the minimum of lij[s)/2 over all the particles 
j for any lij{s) > (see Fig. This leads to 



his)— min Zi,(s)/2= min — — . (13) 

/ij(s)>o s-fi,>o 2s ■ rij 

Substituting into Eq. (|lip . we prove that the Voronoi 
volume is indeed the volume available per particle as cal- 
culated in Eq. (fTU)) : 



min (—i^) ds = Wf. (14) 

s-rii>0 2s ■ ri ' 



Formula (|14p can be rewritten as 
1 



W7°'' = 



ds 



Wtds = (Wn. 



(15) 
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where we define the orientational volume for the i particle 
in the s direction as: 

Wf = Vn min . (16) 

Let us recapitulate and recall the three volumes de- 
fined so far which are interrelated: the Voronoi and the 
available volume which satisfy W™"^ = Wf , and the ori- 
entational volume which satisfies (Wf )s = W™"^- 
quantities are additive, thus they provide the total vol- 
ume of the system: 

i i i 

For isotropic packings, (without the average over s) is 
also additive since the choice of orientation s is arbitrary. 
Thus, we obtain: 

w = Y.^w!)s - (E ^t)s = E ^t- (18) 

i i i 

This property reduces the calculations, since there is no 
need for an orientational average. We define the orienta- 
tional free volume function as: 

W^Wt^Vg, (19) 
and the reduced orientational free volume function 



with its average value, w, over the particles, i, for 
isotropic systems as: 

w = {Wi)i = - — — = ( mm — ) -1. 

^ Vg \\2R s-f,,>o s ■ r,J /. 

(20) 

This average orientational volume function requires only 
averaging over the particles i but not over s. The more 
general form of Eq. ([14]) allows study of anisotropic sys- 
tems, a case left to future work. We notice that the av- 
erage of the orientational volume over the particles (for 
a fixed s) is equal to the average of the Voronoi volume 
over the particles: (W/), = (W,™'')*- 

The relevance of Eq. (fT4|) is that: 

(i) It provides a formula for the Voronoi volume for 
any dimension in terms of particle positions fi. 

(a) It suggests that the Voronoi volume is the volume 
function of the system since it is indeed Wf as calculated 
from an independent point of view. 

(in) It allows for the calculation of macroscopic ob- 
servables via statistical mechanics. 

However, further analytical developments are difficult 
since: 



(i) The volume function is a complicated non-local, not 
pair-wise function of the coordinates. 

(ii) It requires the use of field-theoretical methods. 

(in) It cannot be factorized into a single particle par- 
tition function, implying that there are intrinsic strong 
correlations in the system. Such correlations are implicit 
in the global minimization in Eq. (I14|) which, in prac- 
tice, is restricted to a few coordination shells and defines 
a mesoscopic Voronoi length scale, which will be of use 
below. 

To circumvent the above difficulties, we present a the- 
ory of volume fiuctuations to coarse grain W™'' over this 
mesoscopic length scale. This coarsening reduces the de- 
grees of freedom to one variable, the coordination num- 
ber of the grains, and calculates the average mesoscopic 
volume function Eq. (PO)) amenable to statistical calcu- 
lations. 



IV. THE PROBABILITY DISTRIBUTION OF 
VORONOI VOLUMES 

Next, we develop a statistical theory for the proba- 
bility to find a Voronoi volume in order to calculate the 
mesoscopic volume function by averaging the single grain 
function like in Eq. (|^D|) . For a given grain, i, the cal- 
culation reduces to finding the ball j with min^ Vj / cos 9j 
where the minimization is over all grains (see Fig. [3] for 
notation) . We consider Vj = rij , cos 9j = s ■ fij , and we 
set 2R = 1 for simplicity. While the form Eq. (|20|1 is 
valid for any dimension, the following statistical theory 
differs for each dimension. In what follows we work in 
3d. Results in 2d and large dimension will be presented 
in future papers. 

We consider the minimal particle contributing to the 
Voronoi volume along the s direction at (r, 9) (see Fig. 
[31) with 

c = r/ cos 6*. (21) 

We call this particle, the Voronoi particle. Then the 
quantity to compute is the probability to find all the 
remaining particles in the packing at a distance larger 
than c from particle i along the s direction. That is, we 
compute the inverse cumulative distribution function, de- 
noted P>{c), to find all the grains j with Vj/ cos9j > c, 
and therefore not contributing to the Voronoi volume. 

In terms of Fig. [3^, the minimal Voronoi particle lo- 
cated at c, determining the Voronoi cell, defines an ex- 
cluded volume, represented by the gray volume in the 
figure, where no other particle can be located (otherwise 
they would contribute to the Voronoi volume). We de- 
note this excluded volume as V*{c). Thus, P>(c) repre- 
sents the probability that the remaining particles in the 
packing have Vj/ cos 9j larger than c and therefore are 
outside the gray volume V*{c) in Fig. [3^. 

If we know this inverse cumulative distribution, then 
the mesoscopic free volume function is obtained as the 
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c — r/ cos 0; ^ 



(a) 



V*ic) 




V- NV, 



FIG. 4: A mapping between hard sphere and ideal gas in 
the one dimensional system. A system of total volume V of 
A'^ rods in Id (hard spheres in 3d with volume Vg) can be 
mapped to a system of N points of total size V — NVg by 
simply removing the size of the balls. This mapping allows to 
calculate P> exactly in Id which is shown to be an exponential 
as in Eq. (|26p in the large TV limit. For higher dimensions, 
we cannot just remove the size of the balls NVg since the 
shape of the balls is important as well. This implies that 
the exponential probability is an approximation to the real 
distribution in 3d, as discussed in the text. 



and 




(b) 



mean value of wf over the probability density 



dP>{c). 



FIG. 3: Schematic illustration of the derivation of Pb and 
Pc ■ We plot a 2d case for simplicity, although the calculation 
applies to 3d. (a) Pb: Background term. The considered 
particle (green) is located in the center, the closest particle in 
the s direction is at (r, 0), and the white area is the excluded 
zone r < 2R for the center of any other grain. For a fixed 
c, the grey area with volume V*{c) is the region of the plane 
{r',9') where r'/cosO' < c. If the particle at {r,0) is the 
Voronoi particle defining the boundary of the Voronoi cell in 
the s direction, then no other particle is in the grey zone in 
the figure. The computation of Pb involves the calculation of 
the probability that all the remaining particles in the packing 
are outside the grey excluded volume, V*(c). (b) Pc: Contact 
term. The calculation of this term involves the probability to 
find all the contact particles away from the red area defined 
by the closest contact particle to the direction s. 



ic' - 1) 



l)t^.C: 



dc 



\c^-l)'-^dc 
dc 



(23) 



(c^ - l)dPy 



[C^ - l)dPy 



The integration in Eq. (|23p ranges from 1 to cx3 with 
respect to c since the minimum distance for a ball is for 
r = 1 and = giving c = 1 and the maximum at 
r ^ oo. When changing variables to (iP>, the limits of 
integration c : [l,oo) correspond to the inverse cumula- 
tive distribution function P> : [1,0]. 

Next, we calculate the inverse cumulative distribution 
P>(c). Considering the Voronoi particle at distance c, 
the remaining balls are in the bulk and in contact with 
the ball i. A crucial step in the calculation is to separate 
the distribution into two contributions: a contact term, 
Pc(c), and a bulk term, Pb(c). This separation is im- 
portant because it allows one to obtain the dependence 
of the average volume function on the coordination num- 
ber. The contact term naturally depends on the number 
of contacting particles z, while the bulk term depends 
on the average w. Since P>{c) represents the probability 
that all the balls in the packing except the Voronoi ball 
are outside the grey volume in Fig. [3^, then the geomet- 
rical interpretation of the contact and bulk term is the 
following: 

• Pb{c) represents the probability that the balls in 
the bulk are located outside the grey volume V* (c) 
in Fig. [3^. 

• Pc{c) represents the probability that the contact 
balls are located outside the boundary of the grey 
area marked in red in Fig. [3|3 and denoted S* (c) . 

An assumption of the theory is that both contributions 
are considered independent. Therefore: 



w ^ «), = (c3) - 1, 



(22) 



P>(c) = Pb(c)Pc(c). 



(24) 
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Notice that the meaning of Eq. p4|) is the following: 
we first assume that the Voronoi particle is located at 
c. Then all remaining particles, including contact and 
bulk particles, are outside the grey excluded zone. This 
results in the multiplication of the probabilities as in Eq. 
(El. 



A. Calculation of Pb{c) and -Pc(c) 

In order to calculate the distributions Pb and Pc we 
apply the following assumptions. 

We treat the general case of calculating the probability 
for N particles in a system of volume V to be outside a 
given volume V* when added at random. In 3d, this 
probability is nontrivial since the volume occupied by 
each ball in the packing should be greater than the size 
of the ball. However, this probability can be calculated 
exactly in the case of Id in the large iV limit. In Id, 
the distribution of possible arrangements of hard-rods 
corresponds to the distribution of ideal gas particles by 
removing the volume occupied by the size of the ball as 
we can see in Fig. H) Such a mapping is exact in Id, 
implying an exponential distribution of the free volume. 

The probability to locate one particle at random out- 
side the volume V* in a system of volume V is 

Py{l) = il-V*/V). 

For N independent particles, we obtain: 

p>(iv) = (1-^7^^)^. 

We set V*/V = 1/x and the particle density p = N/V. 
Then 



\pxV* 



(25) 



In the limit of a large number of particles, a; — > oo, we 
obtain a Boltzmann-like exponential distribution for the 
probability of N particles to be outside a volume V*: 



P>{N) cx exp(-py* 



N ■ 



oo. 



(26) 



While the above derivation is exact in Id, the extension 
to higher dimensions is an approximation, since there 
exist additional geometrical constraints. Even if there is 
a void with enough volume to be occupied by a particle 
(the volume of the void is larger or equal than the size of 
the particle), the constraint imposed by the geometrical 
shape of the particle (which does not exist in Id) might 
prevent the void from being occupied. Nevertheless, in 
what follows, we assume the exponential distribution of 
Eq. (p6|) to be valid in 3d as well. 

The background is assumed to be uniform with a mean 
particle density given by: 



p{w) 



N 



N 



V-NV, 



NV„ 



(27) 



where we define the volume fraction = w + 1. There- 
fore, Pb assumes a Boltzmann-like distribution of the 
form analogous to Eq. l[26|) , 



PB{c)^c^p{~p{w)V*ic) 



(28) 



where 

V*{c) = 2n I e{c-r/cos0)dr 

foo /'7r/2 



J e{c~ r/ COS ey sin. edOdr 



2tt I r 

'l "'0 



arccos(r/c) 

2 / sinddedr 



(29) 



= 27r y (1 - r/c)r^dr 
is the volume of the grey area in Fig. [3^. We obtain: 



V*ic)^Vg[ic'~l)~3il~-J 



Therefore, 

Pb (c) = cxp 



(c3 - 1) -3(1 - 1/c) 



w 



(30) 



(31) 



The derivation of the surface term, Pc(c), is analogous 
to that of the volume term. Pc(c) is assumed to have 
the same exponential form, analogous to the background 
form of Pb(c), Eq. (pS)) . despite not having the large 
number approximation leading to Eq. (|26p : 



Pc(c) =exp(-ps(^)5*(c)). 



(32) 



where 



S*{c) = 27r /" 
Jo 



arccos(l/ c) 



sm 



2^(1--), (33) 



is the excluded area marked in red in Fig. [Sb. 

To define ps{z) we first follow the analogy with the 
bulk density Eq. (|27p . p{w) to obtain the particle surface 
density on the sphere with z contacting particles: 



Psiz) 



An — zSo 



(34) 



where S'occ = 27r J^^^ sin 9dd — 27ra is the surface occu- 
pied by a single contact ball, with a = 1 — \/3/2 a small 
value. 

However, we notice that the analogy with Eq. is 
more difficult to justify here since the large number limit 
is lacking for the surface term: the maximum number 
of contacting spheres, the so-called kissing number, is 
12 and for disordered packings is 6 in average. Thus, 



rather than considering Eq. (|32p as the exact form of the 
surface term, we take the exponential form as a simple 
variational Ansatz, where the surface density ps{z) is the 
variational parameter which has to be corrected from Eq. 
(p4)) due to the small number of the balls on the surface. 

A more physical definition of ps{z) is a mean free 
density representing the inverse of the average of S* (c) : 
(5*(c)). The meaning of {S*{c)) for a given number of 
contacting particles is the following: add z contact par- 
ticles at random around a central sphere. The average of 
the solid angles of the gaps left between nearest neigh- 
bor contacting spheres is (5*(c)) (see Fig. [S]). Indeed we 
obtain: 



{S* 



dc 



dS* 



S* {c)ps exp{~psS* (c)) -^c'c 



PS I ^ S*eM-PsS*)dS* 

JQ 

poo 

PS / S* exp{-psS*)dS* 
Jo 



(35) 



= 1/PS 

Then, the surface density is replaced by the following 
definition: 



1 



PS 



{S* 



(36) 



Under these considerations, in order to estimate the 
value of psiz) we first consider a single particle approxi- 
mation. We calculate the mean of S* for a single particle 
z = 1, which gives (S*) = 2tt, since S* ranges uniformly 
from to Att. We note that the value is different from the 
prediction of Eq. in the large number limit. Since 
we expect ps{z) to be proportional to z, under a single 
particle approximation we find: 



psiz) « -, 



single particle approximation. (37) 



Simulations considering many z contacting particles 
suggest that corrections from the single particle approxi- 
mation arc important for larger z. Indeed, a more precise 
value is obtained from simulations: 



Psiz) 



z \/3 



z > 1, 



(38) 



deviating from the single particle approximation of Eq. 
([37]). This numerical calculation consists of adding z ran- 
dom, non-overlapping, equal-size spheres at the surface 
of a ball. For every fix z, the sphere closest to the direc- 
tion s defines the free angle Q* and the surface S* — 26* 
(see Fig. [5]). Thus, the idea is to measure the mean gap, 
{S*), left by z contacting particles along the s direction. 



FIG. 5: The calculation of the surface term Pc involves the 
mean free surface density for a given z (=4 in this example) 
obtained from the angle 6* . Note that we show the 2d case 
for simplicity, but the 6* corresponds to a solid angle in 3d. 



CO 




FIG. 6: A numerical calculation using packing of spheres con- 
firms that 2nps{z) — 2tt/{S*} (dots) slightly deviates from 
the single particle estimation as z (dash line), and is much 
better approximated by {^/3/2)z with particle size correction 
(solid line). 



As seen in Fig. [Bl simulations show that the more pre- 
cise value, Eq. ([38]) is valid rather than the single particle 
approximation, Eq. (|37p . 

This many-particle result can be explained by adding a 
small correction term from the area occupied by one par- 
ticle into Eq. ([57)1 . In the case of many contact particles, 
z > 1, many body constraints imply that the surface, 
(5*), should be corrected by the solid angle extended 
by a single ball, and at the same time reduced by the 
increasing number of contacting particles. Thus, up to 
first order approximation. 



(5*) « {27T + Socc)/z. 



This analysis provides a correction to the surface term 



8 



which can be approximated as 
Ps{z) 



2tt 

z 

'2^ 



>5'oc 



(1 



27r(l + a) 

z 



(39) 



2tt 2 



z > 1, 



It is clear that the above derivation is by no means ex- 
act. It merely interprets the origin of the correction from 
the single particle approximation in order to obtain a 
proper estimation of the surface density. The obtained 
value agrees very well with the computer simulations re- 
sults of Fig. [6] and Eq. ([38)) . and therefore we use it 
to define Pc{c). In comparison with numerical simula- 
tions, the derived form of Pc{c) compares well as we will 
show in Fig. [S] Furthermore, the predicted average vol- 
ume fraction compares well with experiments on random 
close packings. 

From Eqs. ([32)) . ([33]) and ([39)). we obtain the surface 
term: 



Pc(c) = exp [- V3z{l ~ l/c)/2 



(40) 



and the inverse cumulative distribution of volumes takes 
the form: 



P>(c) >= exp 



.i(,3_,_3a-i,)-f,a-l, 

(41) ■ 



the average volume function in those dimensions under 
similar assumptions as used in 3d. 

To solve Eq. we start from the identity: 



Then we find: 
1 = 
dexp 

where a — 



X X 

I — exp( )dx = 1. 

Jo w w 



'-((c^-l)-«(l--)): 
w \ c / 

if(c3-l)-a(l-l) 
w \ c 



(44) 



(45) 



wz 



V3/2. Or 



/•I 1 

= / -(c^ - l)dexp 
Jo w 



a I — 1 dexp 

lo w \ c 



(c3 



a(l 



1 1 



- (c^-l)-a(l--) 
w \ c 



(46) 



The second integration in the right hand side is equal 
to zero only at it; = or it; — > cx), corresponding to two 
trivial solutions at (p — 1 and (p — 0, respectively. The 
only non-trivial solution happens at a = 0, and therefore 
we find 



w(z) 



W{z) = 



2V3 



(47) 



V. THE MESOSCOPIC VOLUME FUNCTION 

Substituting Eq. ((4T|) into Eq. ([23|l . we obtain a self- 
consistent equation to calculate the average free volume, 
w. 



w= /'(c^-l) dexp[-if(c3-l)-3(l-i)V 
Jo I w\ c / 



V3 1 



Since 



f\c'-l)dcxp[-{c'-l)/w], 
Jo 



(42) 



(43) 



then Eq. p2|) can be solved exactly, leading to an an- 
alytical form of the free average volume function, under 
the approximations of the theory. The fact that we can 
solve this equation exactly is a fortuitous event. It is 
worth mentioning that an analogous analysis performed 
in 2d as well as in infinite dimensions does not lead to 
an analytical solution of the self-consistent equation ([^^ 
and therefore only numerical solutions are possible for 



We arrive at a mesoscopic volume function (plotted 
in Fig. [7^) which is more amenable to a statistical me- 
chanics approach for jammed matter. Equation (I47p is 
a coarse-grained "Hamiltonian" or volume function that 
replaces the microscopic Eq. ((T^ to describe the meso- 
scopic states of jammed matter. While Eq. (|14p is dif- 
ficult to treat analytically in statistical mechanics, the 
advantage of the mesoscopic Eq. ([i7|) is that it can easily 
be incorporated into a partition function since it depends 
only on z instead of all degrees of freedom . 

If the system is fully random, so that we can extend 
the assumption of uniformity from the mesoscopic scales 
to the macroscopic scales, we arrive to an equation of 
state relating = w + 1 with z as: 

(48) 



Z + 2V3' 

Thus, Eq. (gS]), plotted in Fig. [Tb, can be inter- 
preted as a equation of state for fully random jammed 
matter. We recall that it has been obtained under the 
approximation of a packing achieving a completely ran- 
dom state; the proper test of such an equation would be 
to add spheres at random and observe their distribution 
and average values. In a second paper of this series @ we 
will show that it corresponds to the equation of state in 
the limit of infinite compactivity when the system is fully 
randomized. Indeed, Eq. (|48p already provides the den- 
sity of random loose packings and random close packings 
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FIG. 7; (a) Mesoscopic volume function of granular mat- 
ter, 'w{z) versus z. (b) Equation of state relating the volume 
fraction and coordination number in the limit of infinite com- 
pactivity. 



when the coordination number is replaced by the isostatic 
limits of 4 and 6, for infinitely rough and frictionless par- 
ticles, respectively. 

So far, our analysis only includes geometrical con- 
straints, but has not made use of the mechanical con- 
straint of jamming. The jamming condition can be in- 
troduced through the condition of isostaticity which ap- 
plies to the mechanical coordination number, Z , for rigid 
spherical grains [3] . For isostatic packings the mechanical 
coordination number is bounded hyZ = d+ l = 4 and 
Z = 2d = 6 and we obtain the minimum and maximum 
volume fraction from Eq. (gHl) as 0rlp =4/(4-1- 2\/3) w 
0.536 and (/)rcp — 6/(6 + 2a/3) « 0.634 as shown in 
Fig. [TJd. We identify these limits as the random close 
packing (RCP) and random loose packing (RLP) frac- 
tions (the maximum and minimum possible volumes of 
random packings of spherical particles). This result is 
fully explored in Jamming II [6| where the phase space 
of jammed configurations is obtained using Edwards sta- 
tistical mechanics of jammed matter based on Eq. (|T7)) . 
We note that the mechanical coordination, Z, is not the 
same as the geometrical one, z (see [6*] for details). For 
now, it suffices to say that we can obtain the two pack- 
ing limits in the case of fully random systems of infinite 
compactivity. 

The experimental studies of on the free volume ver- 
sus coordination number of grains in tomography studies 
of random packings of spheres support Eq. (|48l) (see Fig. 
6 in [|). 



A. Quasiparticles 

Equation ([T7|) should be interpreted as representing 
quasiparticles with free volume w and coordination num- 



ber z. When a grain jams in a packing, it interacts with 
other grains. The role of this interaction is assumed in 
the calculation of the volume function (|T7)) and is im- 
plicit in the coarse-graining procedure explained above. 
Thus, the quasiparticles can be considered as particles in 
a self-consistent field of surrounding jammed matter. In 
the presence of this field, the volume of the quasiparticles 
depends on the surrounding particles, as expressed in Eq. 
(|47)) . The assembly of quasiparticles can be regarded as 
a set of non-interacting particles (when the number of 
elementary excitations is sufficiently low). The jammed 
system can be considered as an ideal gas of quasiparti- 
cles and a single particle partition function can be used to 
evaluate the ensemble. These ideas are exploited in the 
definition of the partition function leading to the solution 
of the phase diagram discussed in Jamming II ^] . 



VI. NUMERICAL TESTS 

In this section, we wish to test some predictions and 
approximations of the theory and propose the necessary 
improvements where needed. The purpose of this calcu- 
lation is, first, to evaluate the predictions of the theory 
regarding the inverse cumulative volume distributions, 
and, second, to test whether the background distribution 
is independent of the contact distribution by comparing 
Pb(c) X Pc(c) with P>(c). 

It is important to note that the predictions of the meso- 
scopic theory refer to quasiparticles as discussed above. 
The quasiparticles of free volume Eq. (|T7)l of fix coordi- 
nation number can be considered as the building blocks 
or elementary units of jammed matter. To properly study 
their properties and test the predictions and approxima- 
tions of the present theory one should in principle isolate 
the behavior of the quasiparticles and study their statis- 
tical properties such as the distribution of volumes and 
their mean value as a function of coordination number. 
Such a study is being carried out and may lead to a more 
precise solutions than the one presented in the present 
paper. 

Nevertheless, in what follows we take an approximate 
numerical route and study the behavior of quasiparticles 
directly from computer generated packings with Molec- 
ular Dynamics. While such packings already contain 
the ensemble average of the quasiparticle, we argue that 
in some limits they could provide statistics for isolated 
quasiparticles, at least approximately. This is due to a 
result that needs to wait until Jamming II [6], where 
we find that for a system of frictionless packings there 
is a unique state of jamming at the mesoscopic level and 
therefore the compactivity and the ensemble average does 
not play a role, at least under the mesoscopic approxima- 
tion developed here. Therefore, below, we use frictionless 
packings in order to test the theory. We note, though, 
that the conclusions of this section remain approxima- 
tive, waiting for a more precise study of the packings of 
fixed coordination number. 
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FIG. 8: Comparison between theory and simulations for the 
inverse cumulative distributions, Pb{c), Pc{c), Pb{c) x -Pc(c) 
and P> (c) for a packing at the frictionless point with z = 6 as 
explained in Jamming II Q . The choice of a frictionless pack- 
ing to test the distribution of volumes is due to the fact that 
these packings are independent of the compactivity as will be 
shown in Jamming II. In general such an equation for the dis- 
tribution of c should be tested by generating local packings 
by randomly adding z balls surrounding a given sphere. 



We prepare a frictionless packing at the jamming tran- 
sition with methods explained in Jamming II. The pack- 
ing consists of 10,000 spherical particles interacting with 
Hertz forces. In the simulation, we pick up a direction s 
randomly, and collect Cb and Cc from the balls, as follows: 



Cb 



s-f ij>o 2s ■ fi 



s-fij>Q 2s ■ fi 



> 2R, 



-,r,,- < 2R. 



(49) 



a.s w = (w^) = (c'^) — 1. While some deviations are 
found in the inverse cumulative distribution, we find that 
the average value of the volumes are well approximated 
by the theory. Indeed, we find that the deviations from 
the theoretical probabilities for u;** > {w'^) appear not 
to contribute significantly towards the average volume 
function. 

For instance, the packing in Fig. [5] has a volume frac- 
tion (f) — 0.64 as measured from the particle positions. 
This value agrees with the average (w*) obtained from 
the prediction of the probability distribution P> (c) . We 
find (c^) = 1.561, then {w") = (c^) - 1 = 0.561 and 
(j) = l/(c^) = l/iiw") + 1) = 0.641 in agreement with the 
volume fraction of the entire packing obtained from the 
position of all the balls, 0.64. 

Thus, the present theory gives a good approximation 
to the average Voronoi volume needed for the meso- 
scopic volume function, even though the full distribution 
presents deviations from the theory. In order to capture 
all the moments of the distribution a more refined the- 
ory is needed. Such a theory will include the corrections 
to the exponential forms of -Pb(c) and Pc{c) and their 
correlations. The main result of the mesoscopic theory, 
being the average Voronoi volume decreasing with the 
number of contacts, is not affected by the assumptions of 
the theory for the full probability distribution. 

The correlations between the contact and bulk term 
are quantified by comparing Pb (c) x Pq (c) with P> (c) in 
Fig. [51 From the figure we see that below and around 
the mean (c) , the full distribution is close to the theoret- 
ical result while deviations appear for larger c. Further 
testing of the existence of correlations between Pb (c) and 
Pc (c) is obtained by calculating the product- moment co- 
efficient of Pearson's correlation as follows. 

The Pearson's coefficient is: 



Thus, for a given direction, we find two minimum val- 
ues of c independently as ct and Cc over all particles in the 
packing. We then collect data for 400 different directions. 
The Cb is only provided by the balls in the background, 
and Cc is only provided by the balls in contact. From 
the probability density, we then calculate the cumulative 
probability to find a ball at position r and 9 such that 
r/ cos^? > c. That is, we calculate the cumulative distri- 
bution of Cb and Cc individually, i.e., PB(cb) and Pc(cc). 

The inverse cumulative distributions are plotted Fig. 
|8] showing that the theory approximately captures the 
trend of these functions but deviations exist as well, es- 
pecially for c values larger than its average. The contact 
term Pq (c) is well approximated by the theory support- 
ing the approximations involved in the calculation of the 
surface density term, Eq. p9p . The background term 
shows deviations for larger c; for smaller c the theory is 
not too far from simulations. From Fig. [5] we notice that 
P>(c) is well reproduced by the theory up to a value of 
10% of its peak value. 

It is important to note that the mesoscopic volume 
function, w, is extracted from the mean value of (w^) 



Sbc 



SbbSc. 



(50) 



where Sbb = S(c6^ - c^). Sec = S(cc^ - c^), and 
Scb = S(cfcCc — CbCc). We find that the Pearson coef- 
ficient — 0.0173 is close to zero, meaning that the 
correlations between Pb{c) and Pc(c) are weak. 

The present numerical results imply that the current 
assumptions of the theory are reasonable. The conclu- 
sions are that while the cumulative distributions present 
deviations from the theory in their tails, the average value 
of the Voronoi volumes are well captured by the approx- 
imations of the theory, therefore providing an accurate 
value for the volume function. 

More importantly, the present approach indicates a 
way to improve the theory to provide more accurate re- 
sults. Our current studies indicate that an exact solution 
of the distribution, P> (c) , may be possible up to the sec- 
ond coordination shell of particles, for a fixed z-ensemble. 
Due to the fact that the range of Voronoi cell is finite, 
it is possible to work out a description for the finite, but 
large, number of degrees of freedom for both disordered 
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and ordered packings through computational hnear pro- 
gramming, in principle. This approach is related to the 
Hales' proof of the Kepler conjecture The present 
theory is a mean-field version in terms of the restricted 
description of the disordered packings, which allows us to 
reduce the dimensionality of the original problem in order 
to write down the analytic form of the volume function 
in reasonable agreement with known values of RCP and 
RLP. The present approximations of the theory are fur- 
ther supported by agreement between the obtained form 
of the volume function and the empirical findings of the 
experiments of @. In a future paper we will present a 
more exact theory of the volume fluctuations capturing 
not only the mean value but also higher moments. 

It is of interest to test the formula Eq. (j47|l with the 
well-known example of the FCC lattice at z = 12 to as- 
sess the approximations of the theory. At this limiting 
number of neighbors the entire class of attainable orienta- 
tional Voronoi cells have volumes in a very narrow range 
around 0.7Vg which is larger than the prediction from Eq. 
(|i7)) . The free volume of the FCC Voronoi cell is 0.35135 
while the mesoscopic volume function for z = 12 gives 
2\/3/12 = 0.2886, below the real value. We explain this 
discrepancy since the current theory is developed under 
the assumption of isotropic packings. Isotropic packings 
are explicitly taken into account in the theory considering 
the orientational Voronoi volume >Vf (along a direction 
s) as a simplification of the full Voronoi volume, W™''. 
Such a simplification is meaningful for isotropic packings 
but fails for anisotropic or ordered packings. Indeed, in 
the case of packings with strong angular correlations, the 
"weak-coupling" hypothesis employed here does not work 
well. The extension of the current theory to anisotropic 
packings, such as the FCC lattice aX z = 12, can be car- 
ried out, but remains outside of the scope of the present 
work. In this case, the full Voronoi volume of Eq. ([M]) 
must be used. 

Eventually a volume function that accounts for dis- 
ordered and ordered packings is very important. This 
volume function could test the existence of phase tran- 



sition between the ordered and disordered phases. The 
theory we propose here is a quasiparticle version in terms 
of the restricted description of the disordered packings, 
which allows us to reduce the dimensionality of the orig- 
inal problem in order to write down the analytic form 
of the volume function. The fact that Eq. (j47|) does 
not predict correctly the volume of FCC but does pre- 
dict correctly the volume of RCP raises the interesting 
possibility that there could be a phase transition at RCP, 
an intriguing possibility which is being explored at the 
moment. 



VII. SUMMARY 

In summary, we present a plausible volume function 
describing the states of jammed matter. The definition 
of Wi is purely geometrical and, therefore, can be ex- 
tended to unjammed systems such as colloids at lower 
concentrations. In its microscopic definition, the volume 
function is shown to be the Voronoi volume per particle, 
and an analytical formula is derived for it. However this 
form is still intractable in a statistical mechanics analy- 
sis. We then develop a mesoscopic version of Wi involving 
coarse graining over a few particles that reduces the de- 
grees of freedom to a single variable, z. This renders the 
problem within the reach of analytical calculations. The 
significance of our results is that they allow the devel- 
opment of a statistical mechanics to predict the observ- 
ables by using the volume function in a Boltzmann-like 
probability distribution of states [3] when the analysis is 
supplemented by the jamming condition 'a la Alexander' 
[?]; which we propose in [Q] to be the isostatic condition 
for rigid particles. 
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